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A recursion operator is an integro-differential operator which maps a generalized symmetry 
of a nonlinear PDE to a new symmetry. Therefore, the existence of a recursion operator 
guarantees that the PDE has infinitely many higher-order symmetries, which is a key feature of 
complete integrability. Completely integrable nonlinear PDEs have a bi-Hamiltonian structure 
and a Lax pair; they can also be solved with the inverse scattering transform and admit soliton 
solutions of any order. 

A straightforward method for the symbolic computation of polynomial recursion operators 
of nonlinear PDEs in (1 + 1) dimensions is presented. Based on conserved densities and 
generalized symmetries, a candidate recursion operator is built from a linear combination of 
scaling invariant terms with undetermined coefficients. The candidate recursion operator is 
substituted into its defining equation and the resulting linear system for the undetermined 
coefficients is solved. 

The method is algorithmic and is implemented in Mathematica. The resulting symbolic 
package PDERecursionDperator .m can be used to test the complete integrability of polynomial 
PDEs that can be written as nonlinear evolution equations. With PDERecurslonOperator .m, 
recursion operators were obtained for several well-known nonlinear PDEs from mathematical 
physics and soliton theory. 
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1. Introduction 



Completely integrable nonlinear partial differential equations (PDEs) have a rich 
mathematical structure and many hidden properties. For example, these PDEs have 
infinitely many conservation laws and generalized symmetries of increasing order. 
They have the Painleve property [3], bi-Hamiltonian (sometimes tri-Hamiltonian) 



structures [2], Lax pairs fl'], Backlund and Darboux transformations [441. 15411. etc. 
Completely integrable PDEs can be solved with the Inverse Scattering Transform 
(1ST) Q, i, [13. Apphcation of the 1ST or Hirota's direct method [sil . [4o| allows 
one to construct explicit soliton solutions of any order. While there are numerous 
definitions of complete integrability, Fokas p^] defines an equation as completely 
integrable if and only if it possesses infinitely many generalized symmetries. A 
recursion operator (also called a formal symmetry or a master symmetry) is a linear 
integro-differential operator which links such symmetries. The recursion operator 
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is thus a key tool for proving the existence of an infinite hierarchy of generahzed 
symmetries [s^] and for computing them sequentially. 

The recursion operator story 5J, |59|] starts with the Korteweg-de Vries (KdV) 
equation, 



ut + Quux + uzx = 0, 



(1) 



which is undeniably the most famous completely integrable PDE. The first few 
generalized symmetries (of infinitely many) for the KdV equation are 



(2) 



U5x- 



Note that generalized symmetries depend on the dependent variables of the system 
as well as the x-derivatives of the dependent variables (in contrast to so-called 
Lie-point or geometric symmetries which only depend on the independent and 
dependent variables of the system). 

The KdV equation is a member of a hierarchy of integrable PDEs, which are 
higher-order symmetries of the KdV itself. For example, the Lax equation [iH ]. 
which is the fifth-order member in the hierarchy, is ut + G^^^ = 0. 

Based on the recursion formula [HH] due to Lenard, in 1977 Olver [s^ ] derived an 
explicit recursion operator for the KdV equation, namely 



n 



Dl + Aul + 2u^D-^. 



(3) 



In ([3]), Dx denotes the total derivative with respect to x, D^^ is its left inverse, 
and I is the identity operator. Total derivatives act on differential functions 54l ]. 



i.e. differentiable functions of independent variables, dependent variables, and their 
derivatives up to an arbitrary but fixed order. 

The recursion operator ([3]) allows one to generate an infinite sequence of local 
generalized symmetries of the KdV equation. Indeed, starting from "seed" or "root" 
G^^\ repeated application of the recursion operator ([3]), 



J 



1,2, 



(4) 



produces the symmetries in (j2]) and infinitely many more. 

Analysis of the form of recursion operators like ([3]) reveals that they can be split 
into a (local) differential part, TZq, and a (non-local) integral part TZi. The differ- 
ential operator TZq involves Dx , , etc., acting on monomials in the dependent 
variables. Barring strange cases 4^, the integral operator TZi only involves 



and can be written as the outer product of generalized symmetries and cosymme- 
tries (or conserved covariants) 
then the Lie derivative 54, 53, 



11,|59[. Furthermore, if 7^ is a recursion operator, 
60[ of 7^ with respect to the evolution equation is 
zero. The latter provides an explicit defining equation for the recursion operator. 

For more information on the historyof cornpletely integrable systems and recur- 
sion operators, see [l|, [H, M, Q [13, Q M, S, H S, Eo, M, M, M, 
Based on studies of formal symmetries and recursion operators, researchers have 
compiled lists of integrable PDEs 



While the computation of the Lie derivative of TZ is fairly straightforward, it is 
computationally intensive and prone to error when done by hand. For example, 
the computation of the recursion operators of the Kaup-Kupershmidt equation 
or the Hirota-Satsuma system (see Section H]) may take weeks to complete by 
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hand and might fill dozens of pages. Using a computer algebra system to carry 
out the lengthy computations is recommended, yet, it is not without challenge 
either; systems like Mathematica and Maple are designed to primarily work with 
commutative algebraic structures and the computation and application of recursion 
operators requires non-commutative operations. 

There is a variety of methods 48 1 to construct recursion operators (or master 

l2Ci. I54II , one first finds a bi-Hamiltonian struc- 



17, 18 



symmetries). As shown in 
ture (with Hamiltonians Oi and Q2) for the given evolution equation and then 
constructs the recursion operator as 7^ = provided G2 is invertible; for 

the KdV equation, @i = + AuDx + 2uxl and 02 = Dx form a Hamiltonian 
pair 43] and TZ = ©iGg ^ yields dSj). The Hamiltonians are cosymplectic operators, 
their inverses are symplectic operators [s^. A complicated example of a recursion 
operator (obtained by composing cosymplectic a symplectic operators of a vector 
derivative Schrodinger equation) can be found in [6ll |. 

A recent approach [49] uses the symbolic method of Gelfand and Dickey [13], and 
applies to non-local and non-evolutionary equations such as the Benjamin-Ono and 
Camassa-Holm equations. 

At the cost of generality, we advocate a direct approach which applies to poly- 
nomial evolution equations. In the spirit of work by Bilge llj, we use the scaling 
invariance of the given PDE to build a polynomial candidate recursion operator as 
a linear combination of scaling homogeneous terms with constant undetermined co- 
efficients. The defining equation for the recursion operator is then used to compute 
the undetermined coefficients. 

The goals of our paper are threefold. We present (i) an algorithmic method in a 
language that appeals to specialists and non-specialists alike, (ii) a symbolic pack- 
age in Mathematica to carry out the tedious computations, (iii) a set of carefully 
selected examples to illustrate the algorithm and the code. 

The theory on which our algorithm is based has been covered extensively in the 
literature [III, [13, IH, [Hi, [13 . Our paper focuses on how things work rather than 
on why they work. 

The package PDERecursionOperator .m [^ is part of our symbolic software col- 
lection for the integrability testing of nonlinear PDEs, including algorithms and 
Mathematica codes for the Painleve test [3, [^ and the computation of conser- 
vation laws [1, i, I1IJ1IJ2I, [H, [H, [H, [il, generalized symmetries [H, ^ and 



recursion operators jlQ . l26l |. As a matter of fact, our package PDERecursionOp- 



erator. m builds on the code InvariantsSymmetries .m [27] for the computation 
of conserved densities and generalized symmetries for nonlinear PDEs. The code 
PDERecursionOperator .m automatically computes polynomial recursion operators 
for polynomial systems of nonlinear PDEs in (1 -|- 1) dimensions, i.e. PDEs in one 
space variable x and time t. At present, the coefficients in the PDEs cannot ex- 
plicitly depend on x and t. Our code can find recursion operators with coefficients 
that explicitly depend on powers of x and t as long as the maximal degree of these 
variables is specified. For example, if the maximal degree is set to 1, then the coef- 
ficients will be at most linear in both x and t. An example of a recursion operator 
that explicitly depends on x and t is given in Section [7.21 For extra versatility, the 
code can be used to test polynomial and rational recursion operators found in the 
literature, computed by hand, or with other software packages. 
Drawing on the analogies with the PDE case, we also developed methods, algo- 



13, 3^, [11] and generalized 



rithms, and software to compute conservation laws [3l|, , 

symmetries [s^ of nonlinear differential-difference equations (DDEs). Although 
the algorithm is well-established ^], a Mathematica package that automatically 
computes recursion operators of nonlinear DDEs is still under development. 
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The paper is organized as follows. In Section [2] we briefly discuss our method 
for computing scaling invariance, conserved densities, and generalized symmetries 
(which are essential pieces for the computation of recursion operators). Our method 
for computing and testing recursion operators is discussed in Section[3l In Sectional 
we illustrate the subtleties of the method using the KdV equation, the Kaup- 
Kupershmidt equation, and the Hirota-Satsuma system of coupled KdV (cKdV) 
equations. The details of computing and testing recursion operators are discussed 
in Section [5l Section [6] compares our software package to other software packages 
for computing recursion operators. In Section [7| we give additional examples to 
demonstrate the capabilities of our software. A discussion of the results and future 
generalizations are given in Section [HI The use of the software package PDERecur- 
sionOperator .m ^] is shown in Appendix lAl 



2. Scaling Invariance, Conservation Laws, and Generalized Symmetries 

Consider a polynomial system of evolution equations in (1 + 1) dimensions, 

Ut{x,t) = F{u{x,t),Ux{x,t),U2x{x,t), . . . ,Umx{x,t)), (5) 

where F has M components Fi, . . . , Fm, u(x, t) has M components ui{x, t), . . . , 
UM{x,t) and Uj^ = d^u/dx^. Henceforth we write F(u) although F (typically) 
depends on u and its x-derivatives up to some fixed order m. In the examples, we 
denote the components of u(x, t) as u, v, . . . , w. If present, any parameters in the 
PDEs are assumed to be nonzero and are denoted by Greek letters. 

Our algorithms are based on scaling (or dilation) invariance, a feature common 
to many nonlinear PDEs. If ([5]) is scaling invariant, then quantities like conserved 
densities, fluxes, generalized symmetries, and recursion operators are also scaling 
invariant 



54( 1 . Indeed, since their defining equation must be satisfied on solutions 
of the PDE, these quantities "inherit" the scaling symmetry of the original PDE. 
Thus, scaling invariance provides an elegant way to construct the form of candidate 
densities, generalized symmetries, and recursion operators. It suffices to make linear 
combinations (with constant undetermined coefficients) of scaling-homogeneous 
terms. Inserting the candidates into their defining equations then leads to a linear 
system for the undetermined coefficients. 



2.1 Scaling Invariance and the Computation of Dilation Symmetries 

Many completely integrable nonlinear PDEs are scaling invariant. PDEs that are 
not scaling invariant can be made so by extendingthe set of dependent variables 



with parameters that scale appropriately, see [2^, l30(] for details. 



For example, the KdV equation ([T|) is invariant under the scaling symmetry 

{t, X, u) {X'^t, X^'^x, X^u), (6) 

where A is an arbitrary parameter. Indeed, upon scaling, a common factor can 
be pulled out. Assigning weights (denoted by W) to the variables based on the 
exponents in A and setting W{Dx) = 1 (or equivalently, W{x) = W{D~^) = — 1) 
gives Wiu) = 2 and W{t) = -3 (or W{Dt) = 3). 

The rank of a monomial is its total weight; in the KdV equation, all three terms 
are rank 5. We say that an equation is uniform in rank if every term in the equation 
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has the same rank. Conversely, requiring uniformity in rank in ([1]) yields 

W{u) + W{Dt) = 2W{u) + W(D^) = W{u) + 3W{D^,). (7) 

Hence, after setting W{D^) = 1, one obtains W{u) = 2W{D^) = 2 and W{Dt) = 
3W{Dx) = 3. So, scaling symmetries can be computed with linear algebra. 

2.2 Computation of Conservation Laws 

The first two conservation laws for the KdV equation are 

Dt{u) + D^{3u^ + U2x) = 0, (8) 
A {u^) + A (4tx3 -ul + 2uuxx) = 0, (9) 

were classically known and correspond to the conservation of mass and momentum 
(for water waves). Whitham found the third conservation law, 

Dt{u^ — ^M^) + Dx{\u^ — Quu^ + 3u'^U2x + ^u^^ — UxUsx) = 0, (10) 

which corresponds to Boussinesq's moment of instability. For ([5]), each conservation 
law has the form 

Dtp{u{x,t),u.x{x,t), ...) + DxJ{\iix,t),Ux{x,t), . . . ) = 0, (11) 

where p is the conserved density and J is the associated flux. 

Algorithms for coinputing conserved densities and generalized symmetries are 
described in 26, 27 , [13, [13, 113, 34, 35]. Our code, PDERecursionOperator .m [§], 



uses these algorithms to compute the densities and generalized symmetries needed 
to construct the non-local part of the operator. For the benefit of the reader, we 
present an abbreviated version of these algorithms. 

The KdV equation ([1]) has conserved densities for any even rank. To find the 
conserved density p of rank R = 6, we consider all the terms of the form 

L'f-'^(")V(x,t), l<i<R/W{u), (12) 

where is the total derivative with respect to x. Hence, since W{u) = 2, we have 

D^u = U4x, dIu^ = 2ul + 2uu2x, dIu^ = u^. (13) 



We then remove divergences and divergence equivalent terms 3j] , and take a linear 
combination (with undetermined coefficients) of the remaining terms as the candi- 
date p. Terms are divergence equivalent if and only if they differ by a divergence, for 
instance uu2x and divergence equivalent because uu2x — = Dx{uUx)- 

Divergences are divergence equivalent to zero, such as U4,x = Dx{u^x)- Thus, the 
candidate p of rank i? = 6 is 

p = ClV? + C2Ux- (14) 



To determine the coefficients q, we require that (jlip holds on the solutions of ([5]). 
In other words, we ffi'st compute Dtp and use ([5]) to remove \it,\itx, etc. 
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For the KdV equation, 

Dtp = —{\d>CiV?Ux + 'iCiU^U^x + 12C2n^ + \2c2UUxU2x + '2C2UxU^x), (15) 

after ut-, utx, etc. have been replaced using ([1]). Then, we require that Dtp is a total 
derivative with respect to x. To do so, for each component Uj of u, we apply the 
Euler operator (variational derivative) to Dtp and set the result identically equal 
to zero [281]. The Euler operator for u is defined as 

K=0 

where m is the highest order needed. In our scalar example there is only one 
component (u = n). Using (jlSp . which is of order m = 4, 

Cu{Dtp) = -18(ci + 2c2)uxU2x = 0. (17) 

To find the undetermined coefficients, we consider all monomials in u and its deriva- 
tives as independent, giving a linear system for Cj. For the example, ci + 2c2 = 0, 
and taking ci = 1 and C2 = —\ gives 

p = u^-\ul, (18) 
which is the conserved density in conservation law (jlOp . 



2.3 Computation of Generalized Symmetries 

A generalized symmetry, G(u), leaves the PDF invariant under the replacement 
u — u + eG within order e [5J]. Hence, G must satisfy the linearized equation 

AG = F'(u)[G], (19) 

where F'(u)[G] is the Frechet derivative of F in the direction of G, 

B (9F 
F'(u)[G] = ^ F(u + 6G)L=o = Y^^DIG)-^, (20) 

i=0 *^ 

where m is the order of F. The KdV equation ([1]) has generalized symmetries for any 
odd rank. To find the generalized symmetry of rank R = 7, we again consider the 
terms in ()12p . This time we do not remove the divergences or divergence equivalent 
terms. The candidate generalized symmetry is then the linear combination of the 
monomials generated by (I12p . For the example, where W{u) = 2, 

D^U = U^x, DIv? = QUxU2x + 2uu3x, DxV? = Su^Ux, (21) 

so the candidate generalized symmetry of rank i? = 7 is 



G = CiU^Ux + C2UxU2x + C3UU3X + C4U5X- 



(22) 
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The undetermined coefficients are then found by computing ([T9|) and using ([5]) to 
remove utjUtx^utxx, etc. Thus, continuing with the example we have 

2(2ci - 3C2)ulu2x + 2(ci - 3C3)UU2^ + 2(ci - 3C3)uUxU3a: + (C2 - 20C4)ul^ 

+ (C2 + C3 - 30C4^)u2xU4x + (C3 - WC4)UxU5x = 0. (23) 

Again, considering ah monomials in u and its derivatives as independent gives a 
linear system for q. From ()23p . 

Ci = 30C4, C2 = 20C4, C3 = 10C4. (24) 

Setting C4 = 1, we find 

G = SOu^Ux + 20UxU2x + lOuusi: + U^x, (25) 
which is the fifth-order symmetry G^^^ in ([2]). 



3. Algorithm for Computing Recursion Operators 

A recursion operator, TZ, is a linear integro-differential operator which links gener- 
alized symmetries |54l |. 

G(i+9) = 7^G(^■^ J = 1,2,3,..., (26) 

where g is the gap and G^-'^ is the j-th generalized symmetry. In many cases, g = 1 
because the generalized symmetries differ by a common rank and, starting from 
G^^^, all higher-order symmetries can indeed be consecutively generated with the 
recursion operator. However, there are exceptions [3] where g' = 2 or 3. Examples 
of the latter are given in Sections 14.21 14.31 17-31 and Appendix |Al Inspection of the 
ranks of generalized symmetries usually provides a hint on how to select the gap. 



If 7^ is a recursion operator for ([5]) , then the Lie derivative [35|, l5J, ISSl] of TZ is 
zero, which leads to the following defining equation: 

^ + 7^'[F(u)] + n o F'(u) - F'(u) o 7^ = 0, (27) 

where o denotes a composition of operators, 7?-'[F(u)] is the Frechet derivative of 
TZ in the direction of F, 

7^'[F(u)] = ^ [DlFiu]) (28) 
and F'(u) is the Frechet derivative operator, i.e. a M x M matrix with entries 



where m is the highest order occurring in the right hand side of ([5|). 
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In the scalar case, F = F and u = n, (|29p simphfies into 



(30) 



Rather than solving (|27|) , we will construct a candidate recursion operator and use 
i\27h to determine the unknown coefficients, as shown in the following two steps. 

Step 1 Generate the candidate recursion operator. 

The rank of the recursion operator is determined by the difference in ranks of the 
generalized symmetries the recursion operator actually connects, 

rank TZij = rank cf - rank G^''^ , i,j = 1,...,M, (31) 

where 7?. is an M x M matrix and G has M components. As before, g is the gap 
and typically g = 1. Yet, there are cases where g( = 2 or 3. 
The recursion operator naturally splits into two pieces jlOl |. 

7^ = 7^o + 7^l, (32) 

where TZq is a (local) differential operator and TZi is a (non-local) integral operator. 
The differential operator TZq is a linear combination of terms 

D^,°u'l^4^---u''^'I, ko,ki,...en, (33) 

where the ki are non-negative integers taken so the monomial has the correct rank 
and the operator has been propagated to the right. For example, 

Dlul = D^{D^uI) = D^{uJ + uD^) = U2J + 2u^Z)^ + uDI, (34) 

which, after multiplying the terms by undetermined coefficients, leads to 

TZq = CiU2xI + C2UxDx + csuDI. (35) 

We will assume that the integral operator TZi is a linear combination of terms 

G«z)-i®/:u(p(^')), i,jen, (36) 

of the correct rank [llljlEil. In ([55]) . ® is the matrix outer product, and C^ip'^-'^) 
is the cosymmetry (Euler operator applied to p^^^). To standardize TZi, propagate 
Dx to the left. For example, by integration by parts, D~^UxDx = Uxl — D~^U2xI- 
As shown in [ll[, the integral operator TZi can also be computed as a linear 
combination of the terms 

G^'^ D-^ ^ iJe'N, (37) 

of the correct rank, where il^^^^ is the covariant (Frechet derivative of p^^^). While 
G^D-i O Cu{p^^^) is strictly non-local, G^'^^D^ ^ ^ ip^^^ contains both differential 
and integral terms. Therefore, it is computationally more efficient to build the 
candidate recursion operator using C^{p^^^) instead of Finally, the local and 
non-local operators are added to obtain a candidate recursion operator (j32p . 
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Step 2 Determine the unknown coefficients. 

To determine the unknown coefficients in the recursion operator, we substitute 
the candidate into the defining equation (I27p . After normalizing the form of the 
terms (propagating the Dx through the expression toward the right), we group the 
terms in like powers of u, Ux,Uxx, ■ ■ ■ Dx,Dx, ■ ■ ■ , and -D^^. Requiring that these 
terms vanish identically gives a linear system for the Cj. Solving this linear system 
and substituting the coefficients into the candidate operator gives the recursion 
operator for ([5]). If Cj = for all i, then either the gap g is incorrect or ([5]) does 
not have a recursion operator. 

While the gap {g) is usually 1, 2 or 3, it is not obvious which value to take for g. In 
Sections 14.21 and 14. 31 we give a couple of examples where g = 2. Starting from G^^), 
the recursion operator then generates the higher-order symmetries G^^^ , G^^^ , . . . . 
Analogously, starting from G^^^ the recursion operator produces G^^^,G(^^, and 
so on. In Section [7.31 we show an example where g = 3. Further details on how to 
select the gap are given in Section [5l 

4. Examples 

4.1 The Korteweg-de Vries Equation 

To illustrate the method, we use the KdV equation ([1]). Reversing the sign of t, 

ut = F{u) = 6uux + U3x (38) 

for scalar u{x,t). From ([2]), the difference in ranks of the generalized symmetries is 

rank G^^^ - rank G^^^ = rank G^^^ - rank G^^^ = 2. (39) 

Therefore, we will assume g = 1, and build a recursion operator with rank IZ = 2. 
Thus, the local operator has the terms and D^ul = ul of rank 2. So, the 
candidate differential operator is 

TZo = ciDl + C2ul. (40) 

Using p^^^ = u and G^^^ = Ux, the non-local operator is 

7^l = C3G(i)D-i/:„(p(i)) = C3UxD-^Cu{u) = csUxD-\ (41) 

where we used Cu given in (fT6|) . Thus, the candidate recursion operator is 

TZ = TZo + ni = ciDl + C2uI + C3UxD-\ (42) 

Note that each term in (j42|) indeed has rank 2. 

Now, we separately compute the pieces needed to evaluate (f27l) . Using (f30l) . we 
readily compute 

F'{u) = dI + 6uDx + 6uxl. (43) 

Since the candidate recursion operator (fl2|) is t-independent, we have dlZ/dt = 0. 
Next, using (|28]l and (|12|) we compute 

n'[F{u)] = {6C2UUx + C2U3x)I + (6C3M^. + 6C3UU2X + C3U4x)Dx^- (44) 
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Using (jl2]) and ([ISj) . we compute 

n o F'{u) = ciDl + (6ci + C2)uDI + (18ci + Oi)u^Dl (45) 

and 

o 7^ = ciDl + (6ci + C2)u£)^ + (6ci + 3c2 + C3)uxDI (46) 

+ 3(2C2U^ + C2U2a; + C3U2x)Dx + {l2C2UUx + Qc^UUx 
+ C2U3X + ScsUSa;)/ + (GcsU^ + 6C3UU2X + C3U4,x)D~^. 

Substituting (f44l) . (l45l) and (f46l) into (f27l) and grouping like terms, we find 

(4ci - C2)UxDI + (6ci - C2 - C3)u2xDx + (2ci - 03)^3^/ = 0. (47) 
So, 2ci = C3 and C2 = 2c3. Taking C3 = 2, gives 

7^ = + 4u/ + 2u^D^\ (48) 



which is indeed the recursion operator ([3]) of the KdV equation 53|] 



4.2 T/ie Kaup-Kupershmidt Equation 

Consider the Kaup-Kupershmidt (KK) equation H, 59], 



ut = F{u) = 20u'^Ux + 25uxU2x + I0uu3x + u^x- (49) 

To find the dilation symmetry for (I49[) . we require that all the terms in (149[) have 
the same rank: 



W{u) + W{Dt) = 3W{u) + W{Dx) = 2W{u) + 3WiDx) 

= 2W{u) +3W{Dx) = Wiu) + 5W{Dx). (50) 

If we set W{Dx) = 1, then Wiu) = 2, W{Dt) = 5 and the rank of ^ is 7. 
Using Invar iantsSymmetr ies .m we compute the conserved densities 

p(i)=^i, p^'^^ = -8u^ + 3ul, (51) 

and the generalized symmetries 

G^^)=Ux^ G^'^'>=F{u) = 20u'^Ux + 25uxU2x + Wuu3x + U5x (52) 

of (|49p . We do not show G^'^^ through G^^^ explicitly due to length. From the 
weights above, the ranks of the first six generalized symmetries are 

rank G^^^ = 3, rank G^^^ = 7, rank G^'^^ = 9, 

(53) 

rank G^^^ = 13, rank G^^^ = 15, rank G^*^) = 19. 
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We assume that rank TZ = 6 and g = 2, since rank G^^) - rank G^^) / rank G^^) - 
rank G^^^ but rank G^^^ — rank G^^^ = rank G^^^ — rank G^^^ = 6. Thus, taking all 
terms of the form Dl.u^ € N) such that rank {Dl.u^) = 6 gives 

+ {cqUUx + C7U3x)Dx + (Cs-U^ + CgU^ + CioUli2x + CuUix)!- (54) 

Using the densities and generalized symmetries above, we compute 

G(^)L>-^£„(p(2)) = UxD-^Cui-Su'^ + 3ul) = -6uxD-^{4u'^ + U2x) (55) 

and 

G(2)D-i£„(p(i)) = F{u)D~'Cu{u) = F{u)D-\ (56) 
Thus, the candidate non-local operator is 

Til = Cl2UxD^^{Au^ + U2x)I + CIS [2Qu^Ux + 25UxU2x + WuUSx + 1452;) Dx^ ■ (57) 

Substituting 'R, = TZo+lZi into ([27|) gives 49 linear equations for Cj. Solving yields 



(58) 



Cl = -^Ci2, C2 = -6C12, C3 = C5 = -18ci2, C4 = -16ci2, 

69 49 35 13 

C6 = ^Ci2, C7 = ^Ci2, C8 = ^Ci2, Cg = ^Ci2, 

Cio = -60ci2, Cii = -41ci2, Ci3 = -C12, 

where C12 is arbitrary. Setting C12 = —2, we obtain 

7^ = I?^ + 12mI?^ + muxDl + (36u2 + 49?X2x) -D^ 

+ 5 (24Mnx. + 7ii3x.) L>^. + (32u^ + 69u^ + S2uu2x + 13m4x) I 

+ 2n^D-i (4^2 + U2x) I + 2F{u)D-\ (59) 



which was computed in 59(] as the composition of the cosymplectic and symplectic 
operators of (|49]) . 

Since g = 2 the symmetries are not generated sequentially via (|4|). Actually, 
G(i) and G(2) in ([521) are the "seeds" (or roots) and one must use (j26p . Indeed, 
from G(i) one obtains G^^) = nG^^\G^^'> = TZG^'^\ and so on. From G^^), upon 
repeated application of TZ, one obtains GW,G(6), etc. Thus, using the recursion 
operator one can generate an infinity of generalized symmetries, confirming that 
(l49l) is completely integrable. 



4.3 The Hirota-Satsuma System 

Consider the system of coupled KdV equations due to Hirota and Satsuma [H, 

Ut = Fi{vl) =?,UUx-2vVx + \u2,x, , ^ 

^ (60) 

Vt = -F2(u) = -2>UVx - V3x, 
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which model shaUow water waves. Solving the equations for the weights, 

iwiu) + WiDt) = 2W{u) + 1 = W{u) + 3 = 2Wiv) + 1, 

[W{v) + W{Dt) = W{u) + W{v) + l = W{v)+3, ^ ' 

yields W{u) = W{v) = 2 and W{Dt) = 3. 

The first few conserved densities and generalized symmetries computed with 
InvariantsSymmetries .m [2^ are 



p(i)=n, p(2)^3^2_2^2^ 



q(i) ^ ( "^A q(2) ^ f Fiiu)\ ^ / 3uux - 2vvx + \u2,x 



(62) 



and are not shown explicitly due to length. Based on the above weights, 
rank p^^^ = 2, rank p(^) = 4, (63) 



and 



3\ . ^('}\ /5 

(64) 



rank G^^^ = ' 3 ) ' ^ank G^^) - , ^ , , 



rank G^^) = Q , rank G^^) = Q . 

We first set 5 = 1, so that rank T^jj = 2, = 1,2. If indeed the generalized 
symmetries were linked consecutively, then 

-Jl — ( '^^^x + C2W/ + C3VI + C5UI + C6f / \ /g^s. 

Using ([62]) . we have 

Ri = ci3G<'>D-l » £„(p(l>) = ci3 f"^) » £.(p™)) (66) 

= =-(":)«"«('°)=«("S' 2). (67) 



Substituting TZ = TZq + TZi into (j27j) . we find ci = • • • = C13 = 0. Thus, the choice 
(7 = 1 appears to be incorrect. Noting that the ranks of the symmetries in (j64p differ 
by 2, we repeat the process with g = 2, so that rank TZij = 4, i, j = 1, 2. Thus, 

^ = fSt' St') + C4iG(^)l)-i /:.(p(2)) + c42G(2)d-i £.(p(i)), (68) 

Vv'^0j21 i'V.0j22/ 

where {TZo)ij, i,j = 1,2, are linear combinations (with different undetermined co- 
efficients) of {D^,uD'^,vD'^,UxDx,VxDx,u'^ ,uv,v^ ,U2x,V2x}- For instance. 



(^0)12 = CllI?^ + (C12M + C13V)DI + {cuUx + Ci5Vx)Dx 

+ (cien^ + + cisv^ + C19U2X + C2of2z)-f- (69) 
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Using (|62|) . the first term of TZi in (j68j) is 



7^S^^ = c4iG(i)D^ 1 ® = C41 (^^^^ D-^ (6w/ -4?;/) 



The second term of TZi in ([68]) is 



= c,,G(')D-^ = C42 (^;[;;]) Z?-^ ® (/ 0) 

Substituting the form of 7^ = 7^o + 7^l = + TZ^^^ + 7^f ^ into ([27]), the Unear 
system for has a non-trivial solution. Solving the linear system, we finally obtain 

with 

(7e)ii = + 8nD2 + l2u^L)^ + 8 (^2u^ + - '^v^^ I 

(7^)l2 = - '^vDI - ^Va:Dx - ^ {Auv + ■U2a.) / - ^UxD'^vI, 
{11)21 = -l^v^Dx - 12v2j + 4vxD-\l - 4 (3™^ + v^^) D-\ 
(7^)22 = -iDt - IQuDl - Su^Dx - y - ^VxD-\l. 

The above recursion operator was computed in [59] as the composition of the 
cosymplectic and symplectic operators of (f5U|) . 

In agreement with g = 2, there are two seeds. Using (|26p and starting from G^^^ 
in (I62p . the recursion operator (|70ll generates the infinite sequence of generalized 
symmetries with odd labels. Starting from G^^^, the recursion operator (j70p gener- 
ates the infinite sequence of generalized symmetries with even labels. The existence 
of a recursion operator confirms that (j60p is completely integrable. 



5. Key Algorithms 

In this section, we present details of the algorithm. To illustrate the key algorithms 
in Sections 15.21 and 15.31 we will use the dispersiveless long wave system [l| , 

Ut = Fiiu) = UVx + UxV, 

n ) ^^^^ 

Vt = F2{u) =Ux + VVcc, 

which is used in applications involving shallow water waves. 



May 15, 2009 10:3 International Journal of Computer Mathematics BaldwinHeremanl- 
JCM2009ArXiv 



14 

5.1 Integra- Differential Operators 

Recursion operators are non-commutative by nature and certain rules must be used 
to simplify expressions involving integro-differential operators. While the multipli- 
cation of differential and integral operators is completely described by 

iy,Di = iy;^^, (72) 

the propagation of a differential operator through an expression is trickier. To 
propagate the differential operator to the right, we use Leibniz' rule 

= E Q^'^D^'^ ™ e N, (73) 

where Q is an expression and Q^''^ is the k-th derivative with respect to x of 
Q. Unlike the finite series for a differential operator, Leibniz' rule for an inverse 
differential operator is 

oo 

D-'Q = - Q'D^' + Q"D^' - • • • = ^,{-1^') D-'-\ (74) 

fc=0 

Therefore, rather than dealing with an infinite series, we only use Leibniz' rule for 
the inverse differential operator when there is a differential operator to the right 
of the inverse operator. In such cases we use 

D-'QD:; = QD^-^ - D-'Q'Dr'. (75) 

Repeated application of (fTSl) yields 

n-1 

D-'QD: = ^(-l)'=g«I)r'=-i + (-l)"D-iQ(") J. (76) 

fe=0 

By using these identities, all the terms are either of the form PD!^ or PD~^QI, 
where P and Q are polynomials in u and its x derivatives. 

5.2 Algorithm for Building the Candidate Recursion Operator 

Step 1 Find the dilation symmetry. 

The dilation symmetry is found by requiring that each equation in ([5]) is uniform in 
rank, i.e. every monomial in that equation has the same rank. If ([5|) is not uniform 
in rank we use a trick. In that case, we multiply those terms that are not uniform 
in rank by auxiliary parameters (a, /?,...) with weights. Once the computations 
are finished we set the auxiliary parameters equal to one. 

Since the linear system for the weights is always underdetermined, we set 
W{Dx) = 1 and this (typically) fixes the values for the remaining weights. 

For the example under consideration, (j7ip . we have the linear system 



W{u) + W{Dt) = W{u) + W{v) + 1 = W{u) + W{v) + 1, 
W{v) + W{Dt) = W{u) + 1 = 2W{v) + 1. 
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Thus, W{v) = ^W{u),W{Dt) = ^W{u) + 1, provided W{a^) = 1. If we select 
W{u) = 2, then the scahng symmetry for (j7ip becomes 

{t, X, u, v) {X~H, X~^x, X^u, Xv). (78) 

In the code PDERecursionOperator .m, the user can set the values of weights with 
WeightRules -> {weight [u] -> 2>. 

Step 2 Determine the rank of the recursion operator. 

Since the gap g cannot be determined a priori, we assume g = 1- Should this choice 
not lead to a result, one could set g = 2 or 3. In the code, the user can set the Gap 
to any positive integer value (see Appendix E])- 

To determine the rank of the recursion operator, we compute the first 5 + 1 
generalized symmetries and then use 

rank n,j = rank cf - rank Gf^ , i, j = 1, 2, (79) 

to determine the rank of TZ. Hence, the rank of the recursion operator TZ can be 
represented (in matrix form) as follows 

rank TZ = ( ^^"^^ '^^^ ^^^^^ ^12\ /ggN 
I rank 7Z21 rank 7^22 / 

In exceptional cases, the rank of the recursion operator might be lower (or higher) 
than computed by (j79p . In the code, the user has some additional control over 
the rank of the recursion operator. For example, in an attempt to find a simpler 
recursion operator, the rank of the recursion operator can be shifted down by one by 
setting RankShif t -> -1. Similarly, to increase the rank of the recursion operator 
one can set RankShif t -> 1 (see Appendix lAl). 
For (|7ip . the first two generalized symmetries and their ranks are 

GC) = (;^) . rank G<.) = (^) , (81) 
Then, using (|79l) and ([80]) . 

rank ^ = (^J 'fj ■ (83) 
Step 3 Generate the (local) differential operator TZq. 

Given the rank of the recursion operator, we take a linear combination of 

D^^u'l'u^' ■ ■ ■ u)^' a^'"+' (5^"+^ • • • , fco, ^1, • • • G N, (84) 



where the ki are taken so the monomial has the correct rank, the operator Dx 
has been propagated to the right, and a, /?, . . . are the weighted parameters from 
Step 1 (if present). 
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For ([7T|) . the (local) differential operator is 

_(ciDx + C2Vl C3DI + C4Ul + C5vDrc + CQV"^! + C7Vxl\ /gcs 

^'-[ csl c,Dx + c,ovI )• 

Step 4 Generate the (non-local) integral operator TZi. 

Since the integral operator involves the outer product of generalized symmetries 
and cosymmetries, we compute the conserved densities up to 

max{rank TZij - rank (G^^))^ + Wiuj) + W{Dx)}, i,j = 1,..., M. (86) 

We add W{uj) in ([86]) because the Euler operator C^j decreases the weight of the 
conserved density by the weight of Uj. In most cases, we take a linear combination 
of the terms 

G»D-i®/:u(p(^'^), i,jGN, (87) 

of the correct rank as the candidate non-local operator. However, there are cases 
in which we must take a linear combination of the monomials in each term of type 
(187p with different coefficients. 
For ()7ip . we only need the cosymmetry of density p^^^ = v, 



C^ipW) = (0 /) . 
Hence, 

G«I)-i0£.(p«)=(U 2^f,y (89) 
Thus, the (non-local) integral operator is 

\0 cuVxD^'^J 

Step 5 Add the local and the non-local operators to form TZ. 
The candidate recursion operator is 

7^ = 7^o + 7^l. (91) 

So, the candidate recursion operator for (fTTI) is 



* Cs/ CgDx + CiQVl + Ci2VxD-^ I ' ^ ) 



5.3 Algorithm for Determining the Unknown Coefficients 



Step 1 Compute the terms in the defining equation (21). 
Step 1.1 Compute 7^t = 

The computation of TZt is easy. Since the candidate recursion operator is t- 
independent one has TZt = 0. 
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Step 1.2 Compute 7^'[F(u)]. 

The Frechet derivative of TZ in the direction of F(u) is given in ()28p . Unhke the 
Frechet derivative (j20p of F(u) in the direction of G (used in the computation of 
generaUzed symmetries), TZ and F(u) in ()28p are operators. 
AppUed to the example (j7ip with 2 components, 



7?/rp^„^l _ (mn^)])ii (^'[F(u)])i2 

'^^*WJ- I (^/[F(U)])21 (7^'[F(U)])22 



(93) 



where 

(7^'[F(u)]),, = (Z),^F(U)) i,i = l,2. (94) 



fe=o 



Exphcitly, 

(7^'[F(u)])n =C2(m, + ot,)/, (95) 
(7^'[F(u)])2l = 0, (96) 



(7^'[F(u)])i2 = +C5 {Ux + VVx) Dx + [ciUVx + (2C6 + Ca,)vUx + 2C61' Wo; 

+ CjVV2x + CyV^ + 07^22.^/ + Cu{2UxVx + UV2x + VU2x)D^^ , (97) 

(7^'[F(u)])22 = cio (n^ + vvx) I + C12 {u2x + vv2x + vl) D~\ (98) 
Step 1.3 Compute F'(u). 

Use formula (|29p to compute F'(u). Continuing with example (|7ip . 

F'(u)=^^^ + "^^ "^^I^^^V (99) 

V -^2: vDx+VxIJ 

Step 1.4 Compose TZ and F'(u). 

The composition of the M x M matrices TZ and F'(u) is an order preserving inner 
product of the two matrices. For example (fTTj) . 



VoF'(u) - A^°F'(u))n (7^oF'(u))l2^ ..... 
V(7eoF'(u))2i (7^oF'(u))22j' 



with 



(7e o F'(u))ii = C3DI + (ci + C5)vDl + {2ciVx + C2V^ + C4U + cef^ + C7Vx)Dx 

+ {ciV2x + C2VVx + CiiU^) /, (101) 
(JZ o F'(u))i2 = C3fL>^ + [clU + 3C37;i. + C^V^) + (201^2, + C2UV + 3C3f2x 
+ C4nt' + 2c5VVx + CqV^ + C'jVVx)Dx + {ciU2x + C2VUx + CsfSa; 
+ CiUVx + C^VV2x + CeV^Wi: + C7V^ + CiiVUx)l, (102) 

(7^ o F'(u))2i = c^Dl + (C8 + Clo) z;!),. + (cg + C12) z;,./, (103) 

(7^ o F'(u))22 = C<dVDl + (cgn + 2c9t>a; + Cio-U^) Dx + (cslix + C9-U2a; 

+ CiqVVx + Ci2VVx)l. (104) 

Similarly, 
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F (u) o 7^ - ^(p,^^^ ^ ^^^^ (p,^^^ ^ J , (105) 



with 



(F'(u) o TZ)u = civDI + [ciVx + C2f^ + csu)Dr, + (2c2f'Ux + cs^x) I, (106) 
(F'(u) o 7^)l2 = c^vDl + (cs^a; + C5V^ + cgu) Dl + {c4,uv + 2c5t;?;a; + cei'^ 

+ C-jVVx + Cgti^; + CioUv)Dx + (c4nfa; + 04^2^ + "iCQV^Vx 
+ Cjvl + CjVV2x + CloUVx + Ciot^n^. + CllVUx + Cl2UVx)I 
+ {cilUxVx + CiiTO2:e + Cl2nU2i. + Cl2UxVx) D^'^ , (107) 

(F'(u) o 7^)2l = cii?2 ^ ^ ^ ^ ^^j^ (108) 

(F'(u) o 7e)22 = C3DI + (C5 + eg) uL*^ _^ (^^^^ _^ ^^.^^ _^ ^^^2 _^ ^^^^ 

+ CgUa; + Ciov'^)Dx + (c4Ua; + 206^^2: + C7V2x + 2cioUfa; 

+ CiiUx + Ci2-Ufa;)/ + (cilli2x + Ci2VV2x + Cuvl) . (109) 

Step 1.5 Sum the terms in the defining equation. 

For ([7T|) , summing the terms in the defining equation ([27|) , we find 



CsD^, + C^vDI + (Ci + C7)VxDx + (C4 - cs) uDx + CqV^D, 



X 



+ (C2 - C8 + Cii) n^J + CiV2xI = 0, (110) 



(Ci - Cg) liZ)^ + 2c3VxDI + (C2 - Cio) UuZ^a; + (C5 - Cg + 2ci) ti^Z^a; 

+ C^VVxDx + SC3V2xDx - (Ci - C7) Uga;/ + (C4 - Ciq - C12) W^:/ 
+ (C2 + 2C6 - Cio) WMxI + C5VV2xI + Cjvll + 03^33;/ 

+ (Cll - C12) UV2xD~^ + (Cii - C12) UxVxD~'^ = 0, (111) 



(Ci - Cg) + (C2 - Cio) + (C2 - C12) fa;/ = 0, (112) 



C3D^ + C5VDI + (C4 - Cg) uDx + (C5 + C7 - Cg) VxDx 

+ CqV^Dx + (C4 - Cs - Cio + Cll) Uxl + 2cqVVxI 

+ (C7 - Cg) V2xl + (Cll - C12) U2xD^^ = 0. (113) 

Step 2 Extract the linear system for the undetermined coefficients. 

Group the terms in hke powers of u, u^^, Ua;^;, D^;, D^, .. . and D~^. Then, 

grouping hke terms and setting their coefficients equal to zero yields a linear system 
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for the undetermined coefficients. For (j7ip . we obtain 

Ci =0, C3 = 0, C5 = 0, C6 = 0, Cy = 0, C4 - Cs = 0, 
Ci - Cg = 0, 2ci + C5 - Cg = 0, C5 + C7 - Cg = 0, C2 - Cio = 0, 
C2 + 2C6 - Clo = 0, C2 - Clo = 0, C4 - C8 - Cio + Cn = 0, C2 - C8 + Cii = 0, 

C4 - Cio - C12 = 0, cii - C12 = 0, C2 - C12 = 0, cii - C12 = 0. 

(114) 

Step 3 Solve the linear system and build the recursion operator. 

Solve the Hnear system and substitute the constants into the candidate recursion 
operator. For ([7T]). we find 

Ci = C3 = C5 = C6 = C7 = Cg = 0, 2C2 = C4 = 2cio = 2cii = 2ci2 = C8, (115) 

so taking cs = 1 gives 

7^=f^f w+¥^ViV (116) 

In [5^ this recursion operator was obtained as the composition of the cosymplectic 
and symplectic operators of (j7ip . 

Starting from G^^\ repeated apphcation of (1116P generates an infinite number of 
generahzed symmetries of (fTl]) . estabhshing its completely integrable. 



6. Other Software Packages 

There has been little work on using computer algebra methods to find and test 



recursion operators. In 1987, Fuchssteiner et al. [23|] wrote PASCAL, Maple, and 
Macsyma codes for testing recursion operators. While these packages could verify 
if a recursion operator is correct, they were unable to either generate the form of 
the operator or test a candidate recursion operator with undetermined coefficients. 
Bilge [11] did substantial work on finding recursion operators interactively with 
REDUCE. 



Sanders and Wang |56l . |59[ wrote Maple and Form codes to aid in the compu- 



tation of recursion operators. Their software was used to compute the symplectic, 



cosymplectic, as well as recursion operators of the 39 PDEs listed in [59(] and [6C 

Recently, Meshkov [i^ implemented a package in Maple for investigating com- 
plete integrability from the geometric perspective. If the zero curvature representa- 
tion of the system is known, then his software package can compute the recursion 
operator. To our knowledge, our package PDERecursionOperator .m [9] is the only 
fully automated software package for computing and testing polynomial recursion 
operators of polynomial evolution equations. 



7. Additional Examples 

7.1 The Nonlinear Schrddinger Equation 

For convenience, we write the standard nonlinear Schrodinger equation (NLS), 

iut + Uxx + '^\u\'^u = 0, (117) 
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as the system of two real equations, 

ut = -{uxx + 2u'^v), 
vt = Vxx + 2uv'^, 



(118) 



where v = u and i has been absorbed in t. 

To determine the weights, we assume W{u) = W{v) so that (jllSp has dilation 
symmetry {t,x,u,v) {\~'^t, X^^x, \u, Xv). Hence, W{u) = W{v) = l,W{Dt) = 2, 
and W{Dx) = l, as usual. The first densities and symmetries are 

pW = UV, p(^) = UxV, (119) 

G(i) = G(2) = (""A . (120) 



Thus, rank TZij = 1, i,j = 1, 2, and the candidate local operator is 

n = ( ^^^"^ ^'^^^ ~^ '^^^^^ ~^ ^'^^^ ~^ '^'^^^^ ^ f 121) 

° \c7Dx + {csU + Cgv)I CiqDx + {cuU + Cuv)! J ' 

The candidate non-local operator is 



Substituting TZ = TZq^TZ\ into (|27[) . solving for the undetermined coefficients, and 
setting ci6 = —2, we find 

Y —IvD^ vl —Ux — 2,vD^ ul J 

Starting from "seed" G^^^ , the generalized symmetries can be constructed sequen- 
tially using ([4]). This establishes the complete integrability of (I118p . 

In [5§] , Wang split (jllTh into an alternate system of two real equations by setting 
u = V + iw. Using the cosymplectic and symplectic operators of that system, she 
obtained a recursion operator which is equivalent to (|123p . 



7.2 The Burgers' Equation 



Consider the Burgers' equation 54 1 



Ut = UUx + Uxx, (124) 

which has the dilation symmetry {t,x,u) — > {X~'^t, X^^x, Xu), or W{Dt) = 2, 
W{u) = 1, with W{Dx) = 1. For (fTM]) . 

p(i) = u, G(^) = ux, = uux + uxx- (125) 

Assuming that g = 1, the candidate recursion operator of rank 7^ = 1 is 

n = ciDx + C2ul + C3G(^)L'-^/:„(p(^)) = ciDx + C2ul + c^UxD-^. (126) 
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Using the defining equation (j27|) . we determine that ci = 2c3 and 02 = 03. Taking 



C3 = 5, gives the recursion operator reported in 53(] 

n = D^ + ^{uI + u^D~^) =D^ + ^D^ {uD~^) . (127) 

As expected, starting from G^^), one computes G*^^) = TZG^^\ G^^^ = TZG^"^^ = 
TZ'^G^^\ etc., confirming that (|124p is completely integrable. 
The Burgers' equation also has the recursion operator |53 |. 



n = m+^{xl + D-^) =tD^ + ^{tu + x)I {tu^ + 1) (128) 

which exphcitly depends on x and t. Using W{t) = —2,W{x) = W{D~^) = — 1 
and W{u) = 2, one can readily verify that each term in (|128p has rank —1. 

To find recursion operators hke (jl28p . which depend explicitly on x and t, we can 
again use scahng symmetries to build TZ. However, one must select the maximum 
degree for x and t. For instance, for degree 1 the coefficients in the recursion 
operator will at most depend on x and t (but not on xt, x^, or which are 
quadratic). To control the highest exponent in x and t, in the code the user can 
set MaxExplicitDependency to any non-negative integer value (see Appendix IXj) . 

With MaxExplicitDependency -> 1, the candidate local operator then is 

Ko = citD^ + {C2X + cztu) I. (129) 
The first symmetries that explicitly depend on x and t (of degree 1) are 

G« = l + t.., (5(^) = i(n + xn.)+tnn. + tn... (130) 

Thus, the candidate non-local operator is 

ill = 04(5(1) L>-i/:„(p(i)) = C4 (tn, + 1) D-\ (131) 

Requiring that R = TZq + TZi satisfies the defining equation ([27|) . next solving 
for the constants ci through C4, and finally setting 04 = ^, yields the recursion 
operator (1128p . Using (jl28p . one can construct an additional infinite sequence of 
generalized symmetries. Furthermore, G*-^-* = iZG^^\ G*-^-* = TZG^"^^ = iZTZG^^\ etc. 
Connections between TZ and TZ and their symmetries are discussed in [13] and [59| . 



7.3 The Drinfel'd-Sokolov- Wilson Equation 
Consider the Drinfel'd-Sokolov- Wilson system [ll. [4l|. 

ut = 3vvx, 

(132) 

vt = 2uv^ + u^v + 2t;3^, 

which has static soliton solutions that interact with moving solitons without de- 
formation. The scaling symmetry for (jl32p is {t,x,u,v) {X~^t, X~^x, X^u, X^v). 
Expressed in weights, W{Dt) = 3,W{D^) = 1, and W{u) = Wiv) = 2. The first 
few conserved densities and generalized symmetries are 
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p(i)=n, p(^)=v^ p(3) = ^^3 _ 2^^2 _ 1^2 ^ ^2^ (^33) 
G« = M, G(^) = f ^f^^^^ V (134) 

\Vx J \UxV + 2uVx + 2t>3a; ^ 

and 

—10u'^Ux + 15v'^Ux + 'i0uVVx —2bUxU2x + ^^VxV2x — ^^UU'ix + 'i^VV^x —2u^x 
IQv^Vx + lf)V^Vx + IQuVUx +AbUxV2x+^^VxU2x + 'i^UV'ix + ^^VU^x + ISt'Sx 

(135) 

Despite the fact that 

rank G^^) = {^^ , rank G^^) = Q , (136) 

we can not take g = \ or 2. Surprisingly, for (1132^ we must set g = 2> and rank 'R-ij = 
6, i,j = 1,2. So, the candidate local operator has elements involving D^. For 
example, 

(7^0)11 = CiD^ + {C2U + C5V) DI + (csUa; + CiQVx) + (caU^ + Cqv'^ + 012^22; 
+Cl3t^2a; + CisUv) + (culiSa; + C15V3X + C20'"'Ux + £21^^; + C25VUX + C26VVx) Dx 
+ {ciV? + C7t;^ + C27VU2X + C28Vf2x + C29UxVx + CgU^ + CiiW^ + CiqU4^x 

+C17V4X + CigUt!^ + C22UU2X + C23UV2X + C2Au'^v) I. (137) 

The candidate non-local operator is 



where 



(^l)ll = - ( \cinU^x + '^CiisUxU2x + \cii<jUU3x + \ci20U^Ux 



5 5 5 \ 2 

--^Ci2lV^Ux - -C122VV3X - T,Ci23VxV2x - -CuaUxD'^v'^ 

b 3 2/3 

2 4 5 

+ ■^Ci25UxD~^U2x + gCi2eUxD~^U^ + -Ci27UVVxD~^ , (139) 
(■^1)12 = -2ci28UxD^^V2x - ^Ci29UxD~'^UV + ScisqWxD'^V , (140) 
(■^1)21 = (ci3lV5x + ^Ci32u'^Vx + ^Ci33VU3x + 70134^^^3; 

\ 9 9 D 

5 35 5 \t-. 1 ^ 12 



3 



35 5 \ 2 

ClSS^f^^Sx + —Ci3eVxU2x + -Ci37UxV2x\ - -Ci3sVxD~^v'^ 



2 4 5 

+ ■^Ci3gVxD~'^U2x + -Cl40VxD;^'^U^ + -CuiUVUxD;^^ , (141) 
(■^1)22 = -2cu2VxD~^V2x + 2ci43V3xD~^V + 2ciuUVxD~'^V 

- ^cu^VxD'^uv + ciAiiVUxD^^v. (142) 
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The terms in (j27p fill 160 pages and grouping like terms results in a system of 508 
linear equations for q. Solving these linear equations and setting ci^q = — 9 gives 
the recursion operator 

with 

{TZ)n =dI + QuDI + l^u^Dl + i^u^ - 2lv'^ + '^U2^DI + fsOnn^. - TSuv^ 



35 \^ 3 2 41 13 69 2 111 

"^T"^"^/ V - 12ui; + —UU2x + —Uix + ^t;'V2a; 

141 o\ / 2 15 2 25 

—Vx)l+ (5. Ux + 5..3. - 15.... - 15..3. + -U.U2. 

- %VxV2x + .5x) ^ + ^.x^^^.2x/ - \uxD~^v^l + UxD~^U^I, (144) 



(7^)l2 = -42?jL>^ - 51.^. - f 48.. + ^.2x) - f 33.... + 60... + ^.3.) 



, 2 15 39 3 \ , 

18.-' + 15.... + 6. . + —uv2x + —vu2x + -t^4. y 

9 

— 27vVxD~^vI — iUxD'^^uvI — -..£'"^.2./, (145) 



(7^)2l = -UvDl - Q7vxDl - (iQuv + ^.2.)^^ - (18... + 53... 

219 \ / 9 o 99 99 27 \ 

+ — .3.j-D. - (^46.... + 2. . + 6. + —uv2x + Y'"4. + -:^vu2xy 

15 35 

15.7J3. + bU^Vx + buVUx + 5..3. + 9.5. + —V^Vx + —VxU2x 

+ Y...2.) + - '^VxD~^v'^I + (146) 

and 

('^)22 = -27L»^ - 54.D^ - 108.. - (27.2 + 33.2 _^ ^^^.^I?^ 



2 

/ 135 \ / 2 27 27 2 

- (^54... + 105... + —""3. j-D. - \ 2^uv + —uu2x + -^Ux 

147 27 201 2\ / \ -1 

+ -—..2. + —■"4. + —r'^x M ~ 9 2... + 2.3. + ... vl 



- 3vxD-\vI - ^..D-i.2.1 (147) 

This recursion operator can also be computed |53] by composing the cosymplectic 
and symplectic operators of (jl32p . Since g = 3 the symmetries are not generated 
via dH). Instead, there are three seeds, G^^\G^^\ and given in (jl34p and (jl35p . 
Using ([26]) . from one obtains G^^^ = TZG^^\ G^^^ = IZG^^^ = TZ'^G^^\ and so 
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on. From G^'^\ upon repeated application of 7^, one gets G^^\G^^\ etc., whereas 
generates etc. Thus, the recursion operator generates a threefold 

infinity of generalized symmetries, confirming that (jl32p is completely integrable. 

This example illustrates the importance of computer algebra software in the 
study of integrability, in particular, for the computation of recursion operators. 
The length of the computations makes it impossible to compute the recursion 
operators for all but the simplest systems by hand. 

8. Conclusions and Future Work 

To our knowledge, no one has ever attempted to fully automate an algorithm for 
finding or testing recursion operators. The commutative nature of computer algebra 
systems makes it a non-trivial task to efficiently implement the non-commutative 
rules needed for working with integro-differential operators. 

Based on our recursion operator algorithm and it implementation in PDERe- 
cursionOperator .m, a large class of nonlinear PDEs can be tested for complete 
integrability in a straightforward manner. Currently our code computes polynomial 
recursion operators for polynomial PDEs (with constant coefficients) which can be 
written in evolution form, = F(u, u^;, . . . , Umx)- 

With the tools developed for finding and testing recursion operators, it would 
be possible to extend the algorithm to find master symmetries as well as cosym- 
plectic, symplectic and conjugate recursion operators. A symplectic operator maps 
(generalized) symmetries into cosymmetries, while a cosymplectic operator maps 
cosymmetries into generalized symmetries. Hence, the recursion operator for a sys- 
tem is the composition of the cosymplectic operator and the symplectic operator of 
the system. A conjugate recursion operator maps conserved densities of lower order 
to conserved densities of higher order. The master symmetry can be used to gen- 
erate an infinite hierarchy of time-dependent generalized symmetries. It would be 
worthwhile to add to our Mathematica code an automated test of the "hereditary" 
condition [2^ for recursion operators. 
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Appendix A. Using the Software Package PDERecursionDperator .m 

The package PDERecursionDperator .m has been tested with Mathematica 4.0 
through 7.0 using more than 30 PDEs. The Backus-Naur form of the main function 
(RecursionOperator) is 

{Main Function) RecuzsionOperaitozKEquations), (Functions), 

(Variables), (Parameters), (Options)] 
(Options) — > Verbose — > (Bool) \ WeightsVerbose — > (Bool) \ 
Gap (Positive Integer) \ 

MaxExplicitDependency — > (Nonnegative Integer) \ 
RankShif t (Integer) \ 
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WeightRules {List of Rules) \ 

WeightedParameters — > {List of Weighted Parameters) 
UnknownCoef f icients — > {Symbol) 

{Bool) True | False 

{List of Rules) {weight[u] — ^(/nteyer), weight[v] — >-(/nte(7er), ...} 

When using a PC, place the packages PDERecursionOperator .m and Invari- 
ant s Symmetries .m in a directory, for example, myDirectory on drive C. Start 

a Mathematica notebook session and execute the commands: 

In[l] := SetDirectory["C:\\myDirectory"] ; (* Specify directory *) 

In[2] := Get ["PDERecursionOperator .m"] (* Read in the package *) 

In [3] := RecursionOperator [ (* Burgers' equation *) 

D[u[x,t] ,t]==u[x,t]*D[u[x,t] ,x]+D[u[x,t] ,{x,2}] , 
u[x,t] , {x,t}] 

Out [3] = 

{{{2CsD^ + Csul + Csu^D-^}}} 

We can find a recursion operator for Burgers' equation which explicitly depends 
on X and t (linearly) by using the option MaxExplicitDependency: 

In [4] := RecursionOperator [ (* Burgers' equation *) 

D[u[x,t] ,t]==u[x,t]*D[u[x,t] ,x]+D[u[x,t] ,{x,2}] , 
u[x,t] , {x,t}, MaxExplicitDependency -> 1] 

Out [4] = 

{{{2C5tD^ + C^ix + tu)I + C5(l + tU:,)D-^}}} 

In [5] := RecursionOperator [ (* Potential mKdV equation *) 

D[u[x,t] ,t]==l/3*D[u[x,t] ,x] ~3+D[u[x,t] ,{x,3}] , 
u[x,t] , {x,t}, 

WeightRules -> {weight [u] -> 1}, Gap -> 2] 

{{{3Ci9-D^ + (Ci + 2Cr9ul)I - 2C^^u^D-^U2:cI]}] 

In this example, we must use the WeightRules option to fix the weights and the 
Gap option to set = 2. 

In [6] := RecursionOperator [ (* Diffusion system *) 

{ D[u[x,t] ,t]==D[u[x,t] ,{x,2}]+v[x,t]~2, 

D[v[x,t] ,t]==D[v[x,t] ,{x, 2}] }, 
{ u[x,t] ,v[x,t] >, {x,t>, 
WeightRules -> {weight [u] -> weight [v]}, 
RankShift -> -1 ] 

{{{C^D^, C2D^ + C^vD-^}, {0, C^D^}}} 



In this system of equations, we again use the option WeightRules to fix the weights. 
We also use the option RankShift to force RecursionOperator to search for re- 
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cursion operator of a lower than expected rank. 

The option Verbose prints out a trace of the calculations, while the option 
WeightsVerbose prints out a trace of the calculation of the scaling symmetry. 
If one or more parameters have a weight, these weights can be specified using 
WeightedParameters. The undetermined constants can be set to any variable us- 
ing the option UnknownCoef f icients (the default is Cj). 
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